#include <stdio.h>
#include <math.h>
   
#define TRUE 1
#define FALSE 0

#define imax 125
#define jmax 485

typedef int logical;

#define clr(xx) xx[ix][iy]=0.0

#define for_all_x for (ix=0;ix<=mx+1;ix++)
#define for_all_y for (iy=0;iy<=my+1;iy++)
#define for_all_points for_all_y for_all_x

#define pointx_ok ((ix>0)&&(ix<mx+1))
#define pointy_ok ((iy>0)&&(iy<my+1))
#define point_ok  ((pointx_ok)&&(pointy_ok))

double u[imax][jmax],
       v[imax][jmax],
       d1[imax][jmax],
       d2[imax][jmax],
       d3[imax][jmax],
       d4[imax][jmax],  /* scratch variable */
       yy[imax][jmax],   /* scratch variable in reset1 and invert */
       s[imax][jmax],
       px[imax][jmax],py[imax][jmax],
       a00[imax][jmax],a01[imax][jmax],a02[imax][jmax],
       a10[imax][jmax],a11[imax][jmax],a12[imax][jmax],
       a20[imax][jmax],a21[imax][jmax],a22[imax][jmax],
       bij[imax][jmax],
       ut[imax][jmax],ut_hld[imax][jmax],
       vt[imax][jmax],vt_hld[imax][jmax],
       d1t[imax][jmax],d1t_hld[imax][jmax],
       d2t[imax][jmax],d2t_hld[imax][jmax],
       d3t[imax][jmax],d3t_hld[imax][jmax],
       st[imax][jmax],st_hld[imax][jmax],pcg_result[imax][jmax],
       x[imax],y[jmax],h,delta_t,delta_t_hld;

char file_name[20];
double xmin,ymin,xmax,ymax,
       dump_stepsize,archive_stepsize,
       t_stop,lambda_zero,laplace_coeff,
       boussinesq_coeff,laplace1_coeff,
       density_spread,viscosity_spread,
       surface_tension_level,
       tension_oil_air,tension_oil_water,
       density_water,density_air,
       density_oil,viscosity_water,
       viscosity_air,viscosity_oil,
       current_t,gravity,delta_t_override,
       distance_ignore;

logical use_axisymmetric,use_archive,project_used;
int fluid_select,mx,my,mx_type,my_type;
static int dump_count=0;
static int number_time_steps=0;
static double plot_stepsize=0.00;

double stream_error;


#define unarchive(xx) fscan_status=fscanf(archive_stream," %le\n",&(xx[ix][iy]))
#define unarchive_d(xx) fscan_status=fscanf(archive_stream," %d\n",&(xx))
#define unarchive_r(xx) fscan_status=fscanf(archive_stream," %le\n",&(xx))

/* read archive data from file_name.arc */
void get_archive()
{
FILE *archive_stream;
char archive_file[25];
int ix,iy,i_close,fscan_status;

  sprintf(archive_file,"%s.arc",file_name);
  printf("***************** reading %s *******************\n",archive_file);
  archive_stream = fopen(archive_file,"r");
  if (archive_stream==NULL)
   printf("cannot find the archive file.\n");

  for_all_points
  if (point_ok)
   {
    unarchive(u); unarchive(v);
    unarchive(d1); unarchive(d2); unarchive(d3);
    unarchive(s);
    unarchive(px); unarchive(py);
    unarchive(ut); unarchive(ut_hld);
    unarchive(vt); unarchive(vt_hld);
    unarchive(d1t); unarchive(d1t_hld);
    unarchive(d2t); unarchive(d2t_hld);
    unarchive(d3t); unarchive(d3t_hld);
    unarchive(st); unarchive(st_hld);
    unarchive(pcg_result);
   }

  unarchive_r(current_t);
  unarchive_d(number_time_steps);
  unarchive_r(delta_t);unarchive_r(delta_t_hld);
  
  i_close=fclose(archive_stream);
}   /* archive_data */



#define dprintf(xx) if (use_debug == TRUE) printf(xx)
#define dprintf2(xx,yy) if (use_debug == TRUE) printf(xx,yy)

FILE *init_parameter_stream,*log_file;
int fscan_status;

logical use_debug;

logical read_logical()
{
logical in_logical;
char in_string[20];

  fscan_status = fscanf(init_parameter_stream,"%s%d\n",&in_string,&in_logical);
  if (in_logical == TRUE)
    {
     dprintf( "TRUE\n" );
    }
  else 
    {
     dprintf( "FALSE\n" );
    }

  return in_logical;
}

int read_int()
{
int in_int;

  fscan_status = fscanf(init_parameter_stream,"%*s%d\n",&in_int);
  dprintf2("%5d\n",in_int);

  return in_int;
}

double read_double()
{
double in_double;

  fscan_status = fscanf(init_parameter_stream,"%*s%lf\n",&in_double);
  dprintf2("%10.4f\n",in_double);

  return in_double;
}

/* read data from eno.dat for initializing the simulation */
void init_parameters()
{
char log_file_name[14],s1[25],s2[25];
logical do_contour,use_compression,always_use_LLF,
  never_use_LLF,front_detect,use_speed,always_use_inversion,
  primitive_mass,primitive_velocity,use_color,use_isothermal,
  use_distance,no_mix,exact_riemann_solver,stokes_flow,
  bias_velocity,bias_density,gravity_fix,ab_flag,area_fix,
  use_curvature;
double ac_constant,alpha_water,ac_detect,
  density_spread_prj,ac_difference_constant,acm_detect;
int r_time,r_space,ac_type,projection_type,convection_type,
    low_space,i_max,j_max;

  printf("initializing parameters.  Opening eno.dat ... \n");
  init_parameter_stream = fopen("eno.dat","r");
  if (init_parameter_stream == NULL)
    printf("Error opening file.\n");

  use_debug = FALSE;
  dprintf("use_debug : ");
  use_debug = read_logical();
  dprintf("do_contour : ");
  do_contour = read_logical();
  dprintf("use compression : ");
  use_compression = read_logical();  
  dprintf("always_use_LLF : ");
  always_use_LLF = read_logical();
  dprintf("never_use_LLF : ");
  never_use_LLF = read_logical();
  dprintf("dump_stepsize : ");
  dump_stepsize = read_double();
  dprintf("CFL number (lambda_zero) : ");
  lambda_zero = read_double();
  dprintf("fluid_select : ");
  fluid_select = read_int();
  dprintf("grid points (x) (mx-1): ");
  mx = read_int()+1;
  dprintf("grid points (y) (my-1): ");
  my = read_int()+1;
  dprintf("time order of accuracy : ");
  r_time = read_int();
  dprintf("spatial order MINUS one : ");
  r_space = read_int();
  dprintf("time cut off : ");
  t_stop = read_double();
  dprintf("ac_constant : ");
  ac_constant = read_double();
  dprintf("ac_difference_constant : ");
  ac_difference_constant = read_double();
  dprintf("alpha_water : ");
  alpha_water = read_double();
  dprintf("density_water : ");
  density_water = read_double();
  dprintf("front detect : ");
  front_detect = read_double();
  dprintf("artifical compression type : ");
  ac_type = read_int();
  dprintf("use_speed : ");
  use_speed = read_logical();
  dprintf("low spatial order MINUS one : ");
  low_space = read_int();
  dprintf("acm_detect : ");
  acm_detect = read_double();
  dprintf("always_use_inversion : ");
  always_use_inversion = read_logical();
  dprintf("primitive_mass : ");
  primitive_mass = read_logical();
  dprintf("primitive_velocity : ");
  primitive_velocity = read_logical();
  dprintf("i_max : ");
  i_max = read_int();
  dprintf("j_max : ");
  j_max = read_int();
  dprintf("use_axisymmetric : ");
  use_axisymmetric = read_logical();

  dprintf("use_color : ");
  use_color = read_logical();
  dprintf("use_isothermal : ");
  use_isothermal = read_logical();
  dprintf("use_distance : ");
  use_distance = read_logical();

  dprintf("no_mix : ");
  no_mix = read_logical();

  dprintf("exact_riemann_solver : ");
  exact_riemann_solver = read_logical();

  dprintf("laplace_coeff : ");
  laplace_coeff = read_double();

  dprintf("boussinesq_coeff : ");
  boussinesq_coeff = read_double();

  dprintf("laplace1_coeff : ");
  laplace1_coeff = read_double();

  dprintf("stokes_flow : ");
  stokes_flow = read_logical();

  dprintf("gravity : ");
  gravity = read_double();

  dprintf("projection_type : ");
  projection_type = read_int();

  dprintf("delta_t_override : ");
  delta_t_override = read_double();

  dprintf("density_spread : ");
  density_spread = read_double();

  dprintf("viscosity_water : ");
  viscosity_water = read_double();

  dprintf("surface_tension_level : ");
  surface_tension_level = read_double();

  dprintf("density_air : ");
  density_air = read_double();

  dprintf("viscosity_air : ");
  viscosity_air = read_double();

  dprintf("viscosity_spread : ");
  viscosity_spread = read_double();

  dprintf("density_spread_prj : ");
  density_spread_prj = read_double();
  if ((density_spread_prj < 0.0) || (density_spread_prj > 5.0))
    printf("error density_spread_prj\n");

  dprintf("convection_type : ");
  convection_type = read_int();

  dprintf("bias_velocity : ");
  bias_velocity = read_logical();

  dprintf("bias_density : ");
  bias_density = read_logical();

  dprintf("gravity_fix : ");
  gravity_fix = read_logical();

  dprintf("ab_flag : ");
  ab_flag = read_logical();

  dprintf("area_fix : ");
  area_fix = read_logical();

  dprintf("use_curvature : ");
  use_curvature = read_logical();

  dprintf("file_name : ");
  fscan_status = fscanf(init_parameter_stream,"%*s%s\n",file_name);
  dprintf2("%s\n",file_name);

  dprintf("use_archive : ");
  use_archive = read_logical();

  dprintf("archive_stepsize : ");
  archive_stepsize = read_double();

  dprintf("density_oil : ");
  density_oil = read_double();

  dprintf("viscosity_oil : ");
  viscosity_oil = read_double();

  dprintf("tension_oil_air : ");
  tension_oil_air = read_double();

  dprintf("tension_oil_water : ");
  tension_oil_water = read_double();

  fscan_status = fclose(init_parameter_stream);

  init_parameter_stream = fopen("eno.dat","r");
  sprintf(log_file_name,"%s.log",file_name);
  log_file = fopen(log_file_name,"w");

  while (feof(init_parameter_stream) == 0)
    {
     fscan_status = fscanf(init_parameter_stream,"%s%s",&s1,&s2);
     fscan_status = fprintf(log_file,"%s %s\n",&s1,&s2);
    }
  fscan_status = fclose(log_file);
  fscan_status = fclose(init_parameter_stream);
}
  
void main()
{
int dsize,tbytes,ix,iy;
double biggest_cfl,cfl_ldiff,force_u,force_v,
       cfl_force_u,cfl_force_v,cfl_laplace,
       tension_u,tension_v,cfl_tension_u,cfl_tension_v,
       test_cfl;

  current_t=0.0;
  number_time_steps=0;
  project_used = FALSE;

  init_parameters();
  dsize=sizeof(double);
  tbytes=dsize*imax*jmax*40;
  printf("allocating tbytes: %10d\n",tbytes);

  for_all_points
   {
    clr(px);clr(py);
    clr(u);clr(v);
    clr(d1);clr(d2);clr(d3);
    clr(s);
    clr(ut);clr(vt);
    clr(d1t);clr(d2t);clr(d3t);
    clr(st);
    clr(ut_hld);clr(vt_hld);
    clr(d1t_hld);clr(d2t_hld);clr(d3t_hld);
    clr(st_hld);
   }
  delta_t_hld=0.0;

  xmin=-2.5;xmax=2.5;
  ymin=-7.5;ymax=2.5;
  h=(xmax-xmin)/mx;

  get_archive();

  printf("current_t,number_time_steps,delta_t,delta_t_hld\n");
  printf("%10.5f%5d%10.5f%10.5f\n",current_t,number_time_steps,
   delta_t,delta_t_hld);
} /* main */
